Self-consistent description of Andreev bound states in Josephson quantum dot devices 
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We develop a general perturbative framework based on a superconducting atomic limit for the de- 
scription of Andreev bound states (ABS) in interacting quantum dots connected to superconducting 
leads. A local effective Hamiltonian for dressed ABS, including both the atomic (or molecular) levels 
and the induced proximity effect on the dot is argued to be a natural starting point. A self-consistent 
expansion in single-particle tunneling events is shown to provide accurate results even in regimes 
where the superconducting gap is smaller than the atomic energies, as demonstrated by a compar- 
ison to recent Numerical Renormalization Group calculations. This simple formulation may have 
bearings for interpreting Andreev spectroscopic experiments in superconducting devices, such as 
STM measurements on carbon nanotubes, or radiative emission in optical quantum dots. 



I. INTRODUCTION 

When a quantum dot is connected to superconduct- 
ing electrodes, the proximity effect drastically modifies 
the dot's electronic structure due to the local formation 
of Cooper pairs. The density of states on the dot thus 
exhibits a gap, so that the formation of discrete sub- 
gap states arisesM These Andreev bound states (ABS) 
play certainly an important role as they may contribute 
a large part of the spectral weight! 2 -' and carry most of the 
Josephson currentPA physical understanding of the ABS 
requires to characterize how these states are connected to 
the atomic (or molecular) levels of the uncoupled quan- 
tum dot, and to describe quantitatively their evolution 
as a function of several parameters, such as gate volt- 
age, Coulomb interaction, tunnel couplings, and super- 
conducting gap. Whereas the ABS have been observed in 
metal-superconductor hybrid structures^, no direct spec- 
troscopy has so far been achieved in quantum dot sys- 
tems. Andreev bound states come in pairs, one state 
above and one below the Fermi level, forming a two level 
system. Consequently, recent interest in the spectroscopy 
of the bound states has also been stimulated by proposals 
to use the latter as a qubitP At present, several routes 
have been suggested, such as STM measurements on car- 
bon nanotubes P microwave cavity coupling, 7 visible light 
emission using a Josephon diode p or noise experiments P 

Experimentally, superconducting quantum dots can 
be realized with carbon nanotubes junctions or semi- 
conducting In As islands. It has been shown that 
quantum dots connected to superconducting electrodes 
can be tuned from a Coulomb blockade regime, to a 
Kondo regime^- * 11 * 12 ! to a weakly interacting Fabry- 
Perot regime by changing local gate voltages The 
Josephson current at zero bias and multiple Andreev re- 
flect ionsatfinite bias voltage have been measured in such 
devices P 3 * 14 * 15 * 16 ^The transition from a 0-junction to a ir- 



junction, namely a reversal in the sign of the Josephson 
current p2l has also been been obser ved when a magnetic 
moment forms on the dot! 18 * 19 * 20 * 21 ! As a possible applica- 
tion of superconducting junctions, nano-SQUID devices 
have also been fabricated.^ 

An exact theoretical description of a quantum dot 
coupled to superconducting leads is only possible when 
the Coulomb interaction is fully neglected. Hence 
the interacting single dot system, as described by the 
Anderson model with superconducting electrodes, has 
been so far analyzed by treating the Coulomb interac- 
tion with v arious a nalytical schemes, such as the mean 
field theoryP ? * 23 ' 24 ' the perturbation expansion in the 
Coulomb interactional or in the tunnel coupling.^ Non- 
perturbative calcula tions, using the Non-Crossing Ap- 
proximation (NCA) J2S122I or the functional Renormal- 
ization Group (fRG)p^ as well as numerical simula- 
tions based on the numerical renormalizat ion g roup 
(NRG) IHMmEQlSimi] or Quantum Monte CarkPEIhave 

also been developped. 

None of the analytical approaches mentioned above is 
able to describe entirely the physics of a quantum dot 
coupled to superconducting leads. Whereas lowest or- 
der perturbation expansions in the tunnel coupling will 
hardly capture the proximity effect induced by the elec- 
trodes^, mean field and weak-interaction approaches will 
miss the Kondo effect. NRG calculations on the other 
hand can capture the physics of such a system over a 
wide range of parameters, but are numerically demand- 
ing and not easily portable to more complex molecular 
systems. More importantly, in the view of describing the 
ABS alone, none of these techniques does provide a sim- 
ple physical picture. Henceforth we will develop in this 
paper a new perturbative approach based on an effective 
local Hamiltonian for dressed ABS, that extends the limit 
of large superconducting gap pro posed previo ushEMZl, 
which was used by many authorsP^E 38 ! 39 ' 4 " ' 41 ' 42 ' This 
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approach will illuminate the nature of the ABS in inter- 
acting quantum dots, which can be generally viewed as 
renormalized superconducting atomic states. This will 
provide as well a simple analytical framework that is ac- 
curate in the most relevant cases, and that may thus be 
useful for interpreting future spectroscopic experiments. 
In particular, calculations provided in the proposals of 
Refs. |S] and [5] only qualitatively apply in realistic situ- 
ations where the gap is comparable or smaller that the 
atomic energies, even when the gap is large compared 
to the hybridization to the electrodes. This interesting 
regime is precisely the one that we want to address in 
the present work. In addition, we note that our formal- 
ism, which incorporates the atomic (or molecular) levels 
from the outset, can easily be extended to describe more 
complex systems, as for instance superconducting dou- 
ble quantum dots or molecule s with mo re complicated 
orbital structure (see e.g. Refs! 43 | 44 | 45 | % 

We organize our paper as follows. In Sec. [TTJ the sys- 
tem is mapped onto an effective local Hamiltonian, sim- 
ilarly to the widely used atomic limit, but including the 
proximity effects due to the superconducting leads. In 
Sec. |III[ the perturbation theory around this limit is 
set up and self-consistent equations for the ABS energy 
renormalizations are derived in order to extend the valid- 
ity of the bare perturbative approach. Sec. |IV| illustrates 
how this expansion can describe ABS in superconduct- 
ing quantum dots over a wide ran ge of parameters, by a 
comparison to available NRG dataP^il 



tron with spin a and wave vector k in the lead i = L. R, 
and n a = S a d a . The leads are assumed to be described 
by standard s-wave BCS Hamiltonians Hi with super- 
conducting gaps A i = A e llfi . The phase difference of 
the latter is noted ip = tpi, — tpn- Furthermore, the 
leads are assumed to have flat and symmetric conduc- 
tion bands, i.e. the kinetic energy e- . measured from the 

Fermi level ranges in [-D, D] and the density of states 
is p = l/(2D). We assume ^-independent and sym- 
metric tunneling amplitudes t between the dot and both 
superconducting leads. The dot has a level energy td 
and Coulomb interaction U . Experimentally, the crucial 
characteristic energy scales, namely Coulomb interaction 
U, total hybridization T = 2irt 2 p and gap A , are typi- 
cally all of the same order of magnitude, 20 4| providing a 
challenge for analytical methods. 

The physics of the quantum dot can be described via 
its Green's function 

G d , d (r) = -(T r 9 d (r)tf d (p)) , (2) 

where the Nambu spinor 

has been introduced. Because we will only be interested 
in stationary equilibrium physics, G d d (r) shall be com- 
puted in the Matsubara frequency formalism. 



II. THEORETICAL FORMULATION 
A. Model 

We focus in this paper on a single-level quantum dot 
coupled to superconducting leads, which is relevant ex- 
perimentally for molecular junctions with large single- 
electron level spacing. A simple Hamiltonian able to 
describe such a system is given by the superconducting 
Anderson model 

H = E H * + Hd + E h t< . C 1 ) 

i=L,R i=L,R 

where 

H t = V e r ct c r . - V f A, ct c f - + h.c) 

k.a k 

H d = ^e d dld a + Un ]ni 

a 

H Ti = EK^ + h.c.) ■ 

k,a 

In the above equations, d a is the annihilation operator of 
an electron with spin a on the dot, c- that of an elec- 



B. Effective local Hamiltonian 

As the above Hamiltonian has no exact solution, some 
approximations must be made. Among the physical in- 
gredients we want to include in a non-perturbative way 
is the local pairing on the dot that is crucial for the evo- 
lution of the Andreev bound states. Furthermore, the 
Coulomb interaction shall be taken into account in an ex- 
act manner in order to describe the atomic states faith- 
fully, and to highlight how these are adiabatically con- 
nected to the ABS. However, the usual development in 
weak tunnel coupling t around the atomic limiP-^ is not 
sufficient to describe the proximity effect at lowest order. 
Therefore, we shall consider in what follows an expansion 
around a superconducting atomic limit. 

Such simple solvable limiting case of the model (JT|) is 
often referred to as the limit of large g ap A — > oo, and 
has been discussed previousl y^ 36 * 38 * 3 -^. Expansions for 
finite A have not however been discussed to our knowl- 
edge, and are the topic of this paper. We emphasize from 
the outset (see equation Q below), that the supercon- 
ducting atomic limit as used normally in the literature 
corresponds to the limit D — > oo (i.e. infinite electronic 
bandwidth), taken before A — > oo. The order of the two 
limits is crucial: if the limit A — > oo was to be taken first, 
the dot would be completely decoupled from the leads 
and the proximity effect would be lost, so that the limit 



3 



of infinite gap would reduce to the usual atomic limit. 
As will be shown now, the superconducting atomic limit 
should rather be interpreted as a low frequency expan- 
sion, i.e. a limit where the gap is much larger than the 
characteristic frequencies of the dot. 

We start off by deriving the Green's function defined 
in Eq. {2} using the equations of motion. Thereby, the 
Coulomb interaction U will at first be omitted for the 
sake of clarity. Note that in the end, U will simply give 
an extra contribution which adds to the effective Hamil- 
tonian. Fourier transformation straightforwardly yields 

G d^d( iuj n) = iw n l - e d a z -t 2 ^2 ^^.^.(iWn)^ . (3) 



In the limit u) n -C A, the Green's function Q becomes 
purely static and reduces to 

E^L,fc l (^«) = 2 ^ arctan (K) (e-*** 6 (T) ■ (5) 

k 

Note that the low frequency limit we consider here yields 
a Green's function that indeed depends on the finite 
bandwidth D, and this shows that the limit A — > 00 shall 
not be taken for the proximity effect to survive. In what 
follows, we will therefore keep both D and A finite. Plug- 
ging Eq. ([5]) into the Green's function G d d (iui n ) leads to 
the same result as would have been obtained with the 
effective local Hamiltonian 

= £e d 4d„ ~ (T v e l ^44 + h.c.) , (6) 

where the local pairing amplitude induced by the leads 
on the dot reads 

T v = T- arctan cos (|) . (7) 

which explicitly depends on the ratio D/A. By an ap- 
propriate gauge transformation for the operators d a , it is 
always possible to choose T v e l ^ L ^ R — as shall be 
done from now on. The complete local effective Hamilto- 
nian is obtained when the Coulomb interaction is taken 
into account again. Defining ^ = e d + ^, the energy level 
of the dot is shifted such that the Hamiltonian clearly ex- 
hibits particle-hole symmetry for £ d = 0: 



In Eq. (J3J , LJ n is a fermionic Matsubara frequency, 
G- - (iuJn) the bare Green's function in Nambu space 

ki,ki K ' 

of electrons with a wave vector k in the lead i, and the 
Pauli matrices a a have been introduced. Transforming 
the sum over wave vectors k into an integral over energies 
yields 



(4) 
I 



#cff = E & d °d„ ~ K\ (d\d\ + h.c.) 




The physical interpretation of this effective local 
Hamiltonian is simple. For finite gap, the quantum dot 
is coupled to both the Cooper pairs and the quasipar- 
ticles in the leads. The Cooper pairs, which lie at the 
Fermi level, are responsible for the proximity effect. The 
quasiparticles give rise to conduction electrons excita- 
tions with energies higher than the gap A. In the limit 
u> n <C A , the quasiparticles are far in energy and the cou- 
pling between them and the dot vanishes, which greatly 
simplifies the physics and makes an exact solution possi- 
ble. Yet, as the dot is still coupled to the Cooper pairs at 
the Fermi level, the proximity effect survives with a local 
pairing term proportional to the hybridization T between 
dot and leads. 



C. Spectrum of the effective local Hamiltonian 

As the Coulomb interaction simply yields an extra en- 
ergy shift of U/2 for both empty and doubly occupied 
dot, the eigenvectors and eigenvalues of the local effec- 
tive Hamiltonian (18]) are readily obtained by a Bogoli- 
ubov transformatiorP, in perfect analogy with solution 
of the BCS Hamiltonian. H c g has thus four eigenstates, 
the singly occupied spin 1/2 states | f) an d | I) with 
energy E9 = E9 = £d, and two BCS-like states given by 



> Gr r . r\iv n ) = 2p arctan — — . 

4^ kt ' kl y^uJn 2 + A 2 J y/uj n 2 + A 2 V Ae 
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FIG. 1: Phase diagram of a simple dot with Coulomb interac- 
tion U, energy level ^ and hybridization F to superconduct- 
ing electrodes in the effective local limit. The transition line 
corresponds to E% = E®_ . 



fore to transitions between states with n and n ± 1 elec- 
trons. Hence, the ABS peaks in the spectral function may 
be interpreted as transitions between the superconduct- 
ing atomic levels of the dot {|tr), |+), |— )}, possibly renor- 
malizcd by single-particle tunneling events neglected in 
-ffeff (to be included in the next section). Furthermore, 
transitions from a spin 1 /2 doublet to a spin singlet neces- 
sarily involve an electron exchange between the dot and 
the superconducting leads. As the states |— ) and |+) 
correspond to the superposition of an empty and doubly 
occupied dot, this electron exchange and the final singlet 
states can be understood within the Andreev reflection 
picture. 

Putting everything together, our effective local Hamil- 
tonian in Eq. (|8| describes the energies of the Andreev 
bound states as transition en ergie s from the spin 1 /2 dou- 
blet to the spin singlet statesP^Sl There are thus four An- 
dreev bound states in the large gap limit for the model 
([lj, with energy ±ao and ±6q which read: 



1+) = U \ u) + v*\q) 

R = -«*|U>+«|0> , (9) 

where |0) denotes the empty dot and | ti) the doubly oc- 
cupied dot. The amplitudes u and v can always be cho- 
sen to be real with u — 1/2 W1 + ^/ J ^ + T v 2 and v — 



1/2 W 1 — £d/y id 2 + r^ 2 - The energies corresponding to 

these BCS-like states are E% = U/2 ± ^ d 2 + T v 2 + £ d . 

As E^_ is always larger than E°_, the effective local 
Hamiltonian has two possible ground states: the low en- 
ergy BCS-like state |— ) or the degenerate spin 1/2 dou- 
blet {| t), | [)}■ In the |— ) state, the energy is minimized 
for ip — 0. Thus, the spin singlet phase corresponds to a 
0-junction (a result well known from the weak coupling 
limilP^'). The transition between the singlet phase and 
the spin 1/2 doublet takes place at £ d 2 + T 2 V = U 2 /4, 
and Fig. [T] shows the corresponding phase diagram for 
variable T v and U. The state adopted by the quan- 
tum dot in the large gap limit therefore results from a 
competition between the local pairing (induced by the 
proximity effect and characterized by the hybridization 
r) and the Coulomb interaction. 

D. Andreev bound states 

As outlined in the introduction, the coupling to super- 
conducting leads induces a gap in the spectral function 
of the dot, inside which discrete Andreev bound states 
can form. The spectral function of the dot shows there- 
fore sharp peaks, which could be measured by STM 6 
or microwave-^bpticaP' experiments as proposed recently. 
These peaks indicate addition energies at which an elec- 
tron may enter (or leave) the dot, and correspond there- 



a = E*-EZ=*L-y/tf + r v ' (10) 

&o = E% - El = | + ^ d 2 + T v 2 . (11) 

The 0/ir transition corresponds to the crossing of the |— ) 
and \<j) states, which occurs for ao = 0. 

III. PERTURBATION EXPANSION AROUND 
THE EFFECTIVE LOCAL HAMILTONIAN 

A. Perturbation theory 

The effective Hamiltonian is not sufficient to obtain 
satisfying results for all regimes of parameters. First, 
-ffeff only describes the 0/7r-junction transition due to the 
competition between a local moment state (stabilized by 
the Coulomb blockade) and a spin singlet (induced by the 
proximity effect). However, if the Coulomb interaction is 
strong (i.e. U T, below the Kondo temperature), 
the local moment can be screened by the Kondo effect, 
which will compete with the superconducting gap for the 
— 7r transition, so that a typical scaling in the ratio of 
the Kondo temperature to the gap A will appear. Also, 
the Josephson current in the 7r-phase identically vanishes 
from -ffeff, as the spin doublet does not disperse with 
the superconducting phase difference, a limitation of the 
large gap limit. On a more quantitative basis, the ex- 
perimental gap A is usually of the order of a few kelvins, 
which is also the typical scale for both T and U in carbon 
nanotube quantum dot devices. 

In order to extend the description of the quantum dot's 
physics, energy corrections shall be calculated with a per- 
turbation theory around the effective Hamiltonian (J8j> . 
Once these corrections have been obtained, physical ob- 
servables like the Josephson current may be computed 
via the free energy F = — iln(Z), with j3 the inverse 
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temperature. Therefore, it is most convenient to work in 
an action based description, which directly yields the par- 
tition function Z. Following Ref. we first integrate 
over the fermions in the leads. Omitting the resulting 
irrelevant constant, the partition function reads 



Z 

Sdot 



V (^ d ,^ d )e~ SdM with (12) 
J2 ^d, n H T .Gl a .(iLJ n )H^ dtn 



+ / drUd^d^d^d^r) 

>0 



(13) 



where we have introduced the Grassmann Nambu spinors 
at Matsubara frequency u> n = (2n + l)ir/[3, 



— H 



1 x 

73 



dj (t) 
d 4 (r) 



denoting the Grassmann fields associated with electrons 
in the dot by d a and d a . 



The perturbation consists of the terms in Eq. ( 12 I that 



are not contained in the action S e g corresponding to the 
effective local Hamiltonian. A simple identification yields 



S, 



off 



£ dr ^rf CT (r)(|- + e d )d a ( T ) - |igd T (r)^(r) - |r p |^(r)d T (r) + Ud^d^d^d^ , (14) 



S P ert = [ dr f dT , Y J *d{T)H T pl % {T-T')H^ d {T , )+ f dr (\T ^{t^t) + \T ^{^{t)) . (15) 
Jo Jo Jo 



Note that S e g contains the local pairing term derived 
in section [II B | The proximity effect is thus treated non- 
perturbatively (just like the Coulomb interaction), which 
is the crucial ingredient of our analytic approach. The 
perturbation S pcrt simply corresponds to the tunnel cou- 
pling between the dot and the electrodes other than the 
lowest order proximity effect. 

The actual corrections are calculated by expanding the 
partition function to the first order in S'port according to 

Z = J V(W,^)e~ s ^ s ^' t 

w Jv(y,V)e- s ^(l~S pclt + ...) (16) 
which we then identify with 

Z = ^e-^+e-^+e-^- , (17) 



where the renormalized superconducting atomic levels 
E a = E a a + 5E a and E± = E% + 8E ± are obtained from: 



e -fiE. w e -PK (i _ p § E a ) (18) 

e -PE ± „ e -PEl {l _ fj5E±) (19) 



Because the Coulomb interaction is taken into account, 
Wick's theorem cannot be used to calculate Z. Instead, 
expectation values are calculated using Lchmann repre- 
sentation. Explicit calculations may be found in the ap- 
pendix. In the zero temperature limit (3 — > oo, the energy 
corrections are 
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SE a = -i 2 ]T 



1 1 2A 

1 1 'i 

E % + (E° - E°) + Er + (E°_ EO) + E, 



008(f) 



E % +{E\-E») %+(£°-£°), 



SE 4 



2A 



Er-(E°-E*) E, 



E,-{E%-E») 



2\T v \uv 



(20) 
(21) 



6E- 



-* 2 £ 

k.a 



1 2A 

1 ? 

2% - (£° _ £0) + £ 



COB(f) 



E,-(E°_-E°)> l2|r - |M " 



(22) 



with the quasiparticle energy E^ = -^/e^ 2 + A 2 



B. Self-consistent renormalization of the energy 



Eqs. (20 1- (22 1 yield the first corrections to the energy 



levels, so that the bound states energies ao and bo are sim- 
ply shifted by 6a = 5E_ - 6E a and 6b = 6E + - 6E a . Ob- 
viously, these expressions are logarithmically divergent 
when the bound states energies ao and bo approach the 
gap edge, and are therefore only valid as long as e.g. 
a > riog[(D + A)/(A - a )]. In the limit of large 
gap A ^> ao, these corrections to a are thus of the or- 
der Tao/A, so that the small dimensionlcss parameter 
is indeed T/A. However, this peculiar logarithmic de- 
pendence of the ABS energy renormalization shows that 
doing a straightforward 1 /A expansion around the effec- 
tive local Hamiltonian will be rapidly uncontrolled, and 
will have a hard time reproducing the logarithmic sin- 
gularities at A close to ao. For this reason, and also 
because the large gap limit becomes trivial for a finite 
electronic bandwidth, as discussed in section |IIB| it was 
indeed more appropriate to single out in the total action 
all terms left over with respect to the local supercon- 



ducting effective Hamiltonian, see equation (15), and do 
perturbation theory around these. 

Because our lowest-order expansion obviously still 
breaks down when the gap becomes comparable to the 
bound state energy, one would naturally seek to resum 
the lcadind logarithmic divergences in equations ( 20 )- 
(22 1. This can be achieved using a self-consistency condi- 



tion inspired by Brillouin-Wigner perturbation theory,^ 
which allows to extend greatly the regime of validity of 
the perturbative scheme. The resulting self-consistent 
equations that we obtain are: 



6a = 



A 



D 



COS 



E 

+ 2\t v \uv 



dc 



(I) 



1 



E-a(A) E + bo 
2 1 



E + a Q 
1 



E - a(A) E + bo E + a 



(23) 



and 
6b 



D 



A 

E' 



cos 



de 



(!) 



l 



l 



E-b(A) E + b 
-2 1 



E + a 
1 



E - b(A) E + b E + ao 



2\T v \uv 



(24) 



with E = %/e 2 + A 2 , and ao, bo have been defined in Eqs. 



@-(|TT]), with a(A) = a + 6a, 6(A) = b + 6b. Note 
that terms like l/(E + ao) have no self-consistency be- 
cause there are no associated divergences. Eq. (23 1 and 



( 24 1 clearly now hold as long as the renormalized ener- 



gies a(A) and 6(A) are not too close to the gap edge, ±A 
respectively. 



IV. RESULTS 
A. Phase diagram 

We start by discussing the 0— n transition line, by com- 
parison to the numerical renormalization group (NRG) 
data by Bauer et alP. Fig. [2] shows the extension to 
smaller gaps A values of the phase diagram obtained 
with unrenormalized local superconducting states for in- 
finite gap (Fig. [IJ. Even though our perturbative ap- 
proach is fairly simple, the results reproduce nicely the 
NRG data of Refs. H and [Ml The analytically obtained 
phase diagram is indeed identical to the NRG data for 
A > r. For smaller A/T, the Kondo effect sets in, but the 
transition lines remain quantitatively correct for £<j near 
±17/2, with increasing deviations from the NRG calcula- 
tions close to the particle-hole symmetric point ^ = at 
large Coulomb interaction U. In this regime, the 0-phase 
possesses a Kondo singlet ground state. As the leads are 
superconductors, the formation of a Kondo resonance in- 
volves the breaking of Cooper pairs. Therefore, the tran- 
sition is now due to the competition between Tk and the 
superconducting gap A, and should occur at ks Tk A. 

Fig. [3] shows a plot of the transition line for £ d — as 
obtained with Eq. ( p3| (solid curve). The vertical, dot- 
ted line depicts the asymptote in the effective local limit. 



7 



U 




— I 1 1 

A / jtT -» » 

A/%T= 1.0 

A/7ir = 0.3 v 

A / 7tr = 0.05 — X- - 



.-■:sr"x xxx x-^- 
doublet 
j i L 








0.2 0.4 



0.6 



FIG. 2: (Color online) Phase diagram of a simple dot with 
Coulomb interaction U, tunnel coupling V to superconducting 
electrodes with gap A for if = and 71T = 0.2 D. The symbols 
indicate NRG data from Ref. [3] and the various lines our 
results. 
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FIG. 3: (Color online) Transition line between a doublet state 
and the BCS-like state at particle hole symmetry £d = (solid 
curve) for <p = and 7rr = 0.2 D. The vertical dotted line 
corresponds to the transition asymptote in the effective local 
limit at A — > oo. The dots indicate NRG data from Ref. 
[2] and the solid line our result. The inset displays the same 
curves on a logarithmic scale. 



The symbols again correspond to NRG da taPThe Kondo 
temperature is given by T K = 0.182 U y^ST /■KUe" nU / sr 
(see for example Ref. HJ. The inset shows on a log- 
log scale that our approach captures an exponential de- 
cay of the transition line with the Coulomb interaction. 
Nonetheless, the suppression of the BCS-like phase ap- 
pears quantitatively stronger than expected: a factor 4 
instead of 8 is found in the exponential factor of Tk- The 
reason for this is that the vertex renormalizations have 
not been taken into account, as discussed in the context 
of U-NGAP. Far away from the particle-hole symmetric 
limit, our results for the Kondo temperature reproduce 
the lowest-order scaling theory for the infinite-[7 Ander- 
son modef^S and are in relatively good agreement with 
NRG data for all A/T values. 



B. Energy renormalizations at particle hole 
symmetry (£ d = 0) 



While Fig. [2] only indicates the transition line between 
the spin 1/2 doublet and the lowest BCS spin singlet, it 
is also instructive to look at the actual renormalization 
of the energy levels while varying the gap A from large 
values to smaller ones beyond the critical point. Fig. [4] 
indicates the renormalized energies of the two Andreev 
bound states (i.e. the difference between the spin 1/2 
doublet and the two spin singlets energies) for different 
hybridizations T. We note that our results are in quanti- 
tative agreement with the NRG calculations of Yoshioka 
and OhashiP^ Several regions need to be distinguished. 
If the gap A is much larger than the bandwidth D, all 
curves collapse at the value U/2 (left hand side of Fig. 
[4 1, since there is no hybridization with both quasiparti- 
cles and Cooper pairs anymore, and one recovers the bare 
atomic levels. When the gap starts to decrease, the prox- 
imity effect simply splits the two Andreev bound states 



according to equations ( 10 1-( 11 ). When the gap becomes 



of the same order as the typical energy scales of the dot 
ao and &01 the superconducting atomic levels start to mix 
with the electrodes, so that the energies renormalizc in a 
non trivial way. One can see that the transition involv- 
ing the highest BCS states ends up touching the gap edge 
for A w U/2, so that half of the ABS are absorbed into 
the continuum above A, as can be seen in Fig. [5j The 
lowest BCS state follows however a downward renormal- 
ization, until the Fermi level is crossed and the ground 
state becomes the 0-state. The difference in behavior be- 
tween the lowest and highest bound states (the former 
being never allowed to leave the superconducting gap) 
can be tracked into equations (23)- (24 1, where level re- 



pulsion effects from the gap edge occur for the low energy 
level |— ) but are canceled for the high energy level |+), 
which is hence allowed to penetrate into the continuum. 
These considerations unveil how the ABS may be adia- 
batically connected to the atomic (or molecular) levels in 
a complicated fashion. 



Again, our simple analytic approach reproduces the 
NRG resultsP over a vast regime of parameters. Yet, 
some deviations are observed in the Kondo regime: we 
find (for the highest hybridization tiY = 0.005 D) that 
the high energy BCS-like state is not absorbed anymore 
into the continuum of states - an artifact of the limits 
of our perturbative approach. Notice also that the en- 
ergy corrections are too important if the gap becomes 
very small, an effect actually due to our underestima- 
tion of the Kondo temperature at particle-symmetry, as 
discussed previously. Finally Fig. [5] shows that, in the 
limit of vanishing gap, our approach is only valid as long 
as a(A) > —A (as has been mentioned in section III B I , 
because the lowest bound state artificially escapes from 
the gap. The expected saturation of a(A) near —A can 
be restored by adding a further self-consistency for terms 
such as \/{E + ao) in (23 1 (not shown here). 
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FIG. 4: (Color online) Renormalization of the Andreev bound 
state energies as a function of Tk /A (the Kondo temperature 
is given in the text). The dashed curves correspond to the 
high energy bound state b(A), the solid curves correspond 
to a(A). All curves have been calculated for U = 0.005 D 
and £d = 0, with several hybridization values tvT/D = 
0.001,0.002,0.005 (from left to right). Symbols are the NRG 
results obtained in Ref. |3T] 



FIG. 6: (Color online) Renormalization of the Andreev bound 
state energies outside particle-hole symmetry. The dotted 
curves correspond to the high energy bound state 6(A), the 
solid curves correspond to <x(A). All curves have been calcu- 
lated for U — 0.5 D and tvT = 0.05 D, with several level shifts 
id/U = 0.3,0.375,0.4,0.425,0.45. Quantitatively similar re- 
sults were obtained by the NRG in Ref. [31] 
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FIG. 5: (Color online) Same data as in Fig. [4] but normalized 
by the gap. 



C. Energy renormalizations outside particle hole 
symmetry (£ d / 0) 



to the NRG data by Yoshioka and OhashiP The more 
particle-hole symmetry is broken, the more the low en- 
ergy bound state moves away from the gap edge, ensur- 
ing even better convergence of our expansion for a given 
value of r. This can be understood given that this bound 
state corresponds to the transition between |— ) and the 
spin 1/2 doublet: outside particle-hole symmetry, the dot 
either seeks to be as empty as possible (for £<j > 0) or 
as occupied as possible (for ^ < 0). Thus, a BCS-like 
wave function will be favored. As a consequence, the 
Kondo effect (that necessitates a singly occupied dot) is 
less favored. This corresponds to a regime where our 
approximation scheme works at best. 

Further understanding can be gained by looking at the 
energies of the Andreev bound states as a function of ^ 
on Fig. |7J We recover the fact that the high energy bound 
states increases in energy by breaking particle-hole sym- 
metry, whereas the low energy bound state has a decreas- 
ing energy. In addition, Fig. [7J shows that the dispersion 
of both ABS weakens for increasing hybridization. In- 
deed, the more the dot is hybridized with the leads, the 
less the Andreev bound state energy is sensitive to the 
bare values of the dot parameters. 



From an experimental point of view, the position of 
the energy level of the quantum dot is the most control- 
lable parameter of the system (by a simple gate voltage) . 
Therefore, it is important to analyze the evolution of the 
Andreev bound states for different values of 

Fig. [6] illustrates how the energies of the bound states 
scale with A for / and can be favorably compared 



D. Superconducting correlations on the dot 

In order to further analyze the evolution of the states 
in the dot as a function of the parameters in the 
model , we investigate now the superconducting corre- 
lations (did]) on the dot. For the effective local Hamil- 
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FIG. 7: (Color online) Evolution of the Andreev bound state 
energies as a function of the dot's energy level for U = 0.005 D 
and A — U. The hybridization takes several values ttT/D = 
0.001,0.002,0.005. 
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FIG. 8: (Color online) Superconducting correlations as a func- 
tion of the gap A (for U = D/200 and £ d = 0). 



tonian, these correlations are zero in the spin doublet 
phase. In the BCS-like phase, the correlations are max- 
imal if the two states |0) and | ||) are equivalent, i.e. 
at particle hole symmetry. If the dot level is far from 
£d = 0, the wave function will be predominantly |0) (if 
£d is positive) or | fj.) (if ^ is negative). This obviously 
kills the superconducting correlations. 

As the gap decreases from infinity, the (formerly) 
singly occupied state will start having a BCS-like ad- 
mixture and therefore a non zero superconducting corre- 
lation. In contrast, the mixing will result in a decreased 
correlation in the BCS-like phase. Nevertheless, if the 
gap tends to zero, one would expect the correlations to 
vanish as well. This is indeed what Fig. [8] shows. For 
large gaps, the dot is in the spin 1/2 phase; the corre- 
lations are small, but increase as the states mix. The 
transition to the BCS-like phase results in a discontinu- 
ous jump in the correlations, before they finally vanish 
for very small gaps. It can thus be concluded that the 
correlations should be normalized by the gap if one is 
interested in measuring only the mixing effect. Finally, 
the two different curves show how hybridization stabi- 
lizes the BCS-like state with respect to the spin doublet 
via the — n transition. 

As the Coulomb interaction tries to prevent the for- 
mation of a Cooper pair wave function, the transition 
between the BCS-like phase and the spin doublet can 
also be achieved if the Coulomb interaction is tuned, as 
shown in Fig. [9] The effect of the mixing is clearly visible 
by an increase of the correlation (djrfj) (now normalized 
by the gap) while U is lowered. We also find that the 
correlations relative to the gap decrease for higher gaps, 
which is a simple saturation effect (the highest possible 
correlations are (djdj) = 0.5). Furthermore, our results 
are quantitatively precise if the gap if larger than the hy- 
bridization r for all values of U, while relatively small 
deviations appear for A < 7rr, as shown by the compar- 
ison to the NRG data from Ref. 2 and to second-order 
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FIG. 9: (Color online) Superconducting correlations as a func- 
tion of the Coulomb interaction U, for A/(tiT) = 0.3, 1.0 and 
at particle-hole symmetry £d = 0. Solid lines are the results 
of our self-consistent equations, diamonds correspond to NRG 
data from Ref. O and crosses are perturbation theory in U 
at second order. 



perturbation theory in U (valid in the singlet ph ase o nly, 
providing accurate results for U < 27rr roughly) 25 * 51 !. 

Finally, we analyze how the correlations evolve outside 
particle hole symmetry. As mentioned above, one expects 
the correlations to decrease because the dot evolve from 
a superconducting atomic limit toward a usual atomic 
limit (i.e. from the states |±) toward the states |0) and 

t|))- On the other hand, there will be a transition 
from the spin doublet to the singlet phase and therefore 
a mixing effect. Fig. [TO] shows the competition between 
the mixing effect (that increases the correlations outside 
particle hole symmetry) and the evolution toward the 
normal atomic limit (that lowers the correlations) if £d is 
increased. The effect of the Coulomb interaction is once 
more found to favor the single occupancy. 
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particle hole symmetry (for 7rr = 0.2 D, U — 6F and A = 
0.LD). 



FIG. 11: (Color online) Josephson current through the bound 
states for U = 3T, A = 0.1 D. 



E. Josephson current 

We now turn to the Josephson current through the 
quantum dot. The latter is given by J = 2e^ (where F 
is the free energy). At zero temperature, the free energy 
is the same than the level energies, so that the Josephson 
current can readily be obtained once the renormalized 
energy levels have been calculated. 

Nevertheless, our analytical approach only describes 
the effective local limit atomic states, and we can there- 
fore only determine the current through the Andreev 
bound states. Yet, it is known that the Josephson current 
also contains a contribution of the continuum of states P 
The latter can be of the same order and opposite sign as 
the bound state contribution. Furthermore, Bauer at alP 
have shown that the spectral weight of the bound states 
may vary importantly as a function of the different pa- 
rameters (like the Coulomb interaction U), especially in 
the spin doublet phase. As we exclusively investigate the 
effective local limit states, we do not keep track of this ef- 
fect as well. Therefore, the Josephson currents obtained 
in our approach will only provide a rather rough and 
qualitative idea of the actual total Josephson current. 

Fig. [TT] shows the Josephson current calculated as the 
phase derivative of the ground state energy E x , J = 
for different values of One notices two regimes: If 
the phase is close to ip = 0, the system will be in the 
BCS-like state. As there is no magnetic moment in this 
phase, the ground state corresponds to a 0-junction (i.e. 
phase difference ip = 0). If ip increases, the energy of the 
BCS-like state increases (as can be understood in the 

effective local limit, where E- = U/2 — J £ d 2 + T v 2 ). 
When the BCS-like state crosses with the spin doublet, 
the ground state changes and the dot becomes singly oc- 
cupied. This magnetic moment leads to a discontinuous 
jump in the Josephson current and the formation of a 
7T-j unction. Again, we notice that the spin doublet is 
stabilized in the particle hole symmetric case. 



V. CONCLUSION 

In this section we summarize our main results. First, it 
has been shown how the Hamiltonian of a quantum dot 
coupled to superconducting leads can be mapped onto 
an effective local model if the superconducting gap A 
is much bigger than the characteristic energy scales of 
the dot. This limit can be quite generally regarded as a 
low frequency expansion of the Green's function of the 
dot rather than the limit A — > oo used in the literature. 
This enabled us to extend the effective local Hamiltonian 
to leads with a finite electronic bandwidth. 

We have then set up a perturbation theory around this 
local effective Hamiltonian and established self-consistent 
equations for the energy renormalizations of the Andreev 
bound states. We have derived those equations based on 
the fact that the latter correspond to transitions between 
different states of the local effective Hamiltonian. 

In a last section, we used our formalism to calculate 
physical quantities such as the Andreev bound state en- 
ergies or superconducting correlations, and understood 
how these evolve as a function of gate voltage, hybridiza- 
tion, Coulomb interaction and superconducting gap am- 
plitude. It has been shown that our simple approach 
agrees well with NRG data in a vast range of parameters, 
with the main limitation that the Kondo temperature is 
not quantitatively described near particle-hole symmetry. 
However, most experimentally interesting regimes should 
be described correctly by the simple equations we have 
derived. 

The simplicity and portability constitute the main ad- 
vantages of our approach, if one is interested in the An- 
dreev bound states only, compared to extended numerical 
simulations. As the perturbative description is analyti- 
cal and based on atomic-like levels, it should in principle 
be able to describe more complex systems like multiple 
quantum dots or molecules with several orbitals coupled 
to superconducting environments, and be readily appli- 
cable to describe future spectrosopic measurements. Ex- 
tensions of our formalism to the computation of the tun- 
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neling current at realistic gap values in three-terminal 
geometries^ relevant for STM experiments should cer- 
tainly deserve further scrutiny. 
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calculations are performed in the operator formalism. It 
is very useful to note that the product of two fermionic 
(or bosonic) Greens functions G a (r) and Gf,(r) obeys 



jyr^dT'G a (T-T')G b (T-r') = (3 J Q P drG a (r)G b (r) 
(as can be shown using Fourier transformation). The 
partition function's perturbation expansion is 



APPENDIX: DERIVATION OF THE ENERGY 
CORRECTIONS 



The partition function is derived starting from the ac- 



tion's perturbation expansion in section III The actual 



Z = 



Z - Zot 2 ^ J dr (Gl^ n (r)(T T d\(r)d^Q)) 

k,i 

- G k ;12 (-)M(-)4(°)>o - GL ;21 MW(^t(»)>o + G° &&;22 (r)(T r d i (r)4(0)) 
-2p\T v \ (<T T 4(0)4(0))o+ <T T d i (r)d T (0)) 



(A.l) 



In the above equation, G°- - is the Fourier transformed 

kiki-ij 



and the subscript 



Nambu matrix clement G°- - (ioj n ) 

ki 7 ki 

indicates that the expectation values are statistical av- 
erages calculated in the effective local limit. The leads' 
Green's functions are: 



with E? = \/er 2 + A 2 . Furthermore, G°- - (r) 

K V k ' kiki;21 y ' 

G liki;l2^* aild ^fc G L&;2 2 ( r ) = S£ G Lk;1i( T )- 



T^OK 



k 

A 
A 

2% 



-\r\E B 



2 hfl \E) (E )\,^ S ° ne cann0 ^ a PP^ Wick's theorem because of the 
U r l ) n F\ jtj ^Coulomb interaction, the dot's Green's functions are cal- 
culated using Lchmann representation, which yields (for 
r > 0) 



e -\T\E % „ e -{0-\r\)E % 
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(T T 4(r)d T (0)) = ^ (« 2 (e-^-^e-^OJ-^) + e -^ e -n03-^)) 

(T T d\(r)dUO)) = ^UV ( e -^r e -EO_ { 0-r) _ p -E»_r e -E^-r) 

_ e -E°r e -E° + (0-r) + e -E° + Tg-^-r)) (A3) 

(T T rf i (r)d T (0)) o = ^™( e -^ e - B -(^)- e - B -- e - £ ?(^) 

- e - s ?-e- E +^-^ + e-^Ve-^?^-)) (A.4) 
(T T d i (r)4(0)) = ^(n 2 ( e - B +- e - B ?^) + e-^e- B -^)) 



Using u 2 + u 2 = 1, the partition function becomes: 



z = z Q + pt 2 J2 
1 

X 



+ 



E-{E%-E») 
1 

E+(E° + -E°) 
2A 



1 



25 



-P(e+e; 



2) ) 



E + (E°_ -E%) 



E 



-uv 



E + (E% -E%) 
1 



(e~ 0E ° - e-PiE+El)^ 



E-{E%-E% 
+2(3\T v \uv (e-^+-e-^-) 



/3S^ _ t - HE ■ f; 
e -l3E 



£ + (££ - £0) 
E-(E°_- Ej) ( 



e -/3B„ _ e -/3(B+B 



e -/3S° _ e -/3(£+B° 



(A.6) 



As = i/eg 2 + A 2 > 0, terms with an e /3E * are ex- 
ponentially suppressed for T — ► K and can be omitted. 
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